Slow decay of concentration variance due to no-slip walls in chaotic mixing 
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Chaotic mixing in a closed vessel is studied experimentally and numerically in different 2-D flow 
configurations. For a purely hyperbolic phase space, it is well-known that concentration fluctuations 
converge to an eigenmode of the advection-diffusion operator and decay exponentially with time. We 
illustrate how the unstable manifold of hyperbolic periodic points dominates the resulting persistent 
pattern. We show for different physical viscous flows that, in the case of a fully chaotic Poincare 
section, parabolic periodic points at the walls lead to slower (algebraic) decay. A persistent pattern, 
the backbone of which is the unstable manifold of parabolic points, can be observed. However, 
slow stretching at the wall forbids the rapid propagation of stretched filaments throughout the 
whole domain, and hence delays the formation of an eigenmode until it is no longer experimentally 
observable. Inspired by the baker's map, we introduce a 1-D model with a parabolic point that 
gives a good account of the slow decay observed in experiments. We derive a universal decay law 
for such systems parametrized by the rate at which a particle approaches the no-slip wall. 



PACS numbers: 47.52. +j, 05.45.-a 
I. INTRODUCTION 

Many industrial applications involve the mixing of vis- 
cous fluids. Fields as diverse as chemical engineering, the 
pharmaceutical and cosmetics industries, and food pro- 
cessing depend on the stirring of initially heterogeneous 
substances to obtain a product with a sufficient degree of 
homogeneity. Viscous, confined or fragile fluids are best 
mixed by non-turbulent flows, which tend to be less ef- 
fective at mixing than their turbulent counterpart. How- 
ever, some laminar flows exhibit chaotic advection, mean- 
ing that they have chaotic Lagrangian trajectories [1, 2], 
allowing them to rival turbulent flows in their ability to 
mix. The framework of chaotic advection and dynamical 
systems provides a useful characterization of mixers that 
relies on the nature of the phase space, or the stretching 
statistics of Lagrangian trajectories. However, an essen- 
tial issue for the mixing of a diffusive passive scalar is 
to predict the rate at which the scalar concentration is 
homogenized by a given stirring protocol. 

Various approaches including an eigenmode analy- 
sis [3-12], a large-deviation description of the stretching 
distribution [13, 14], and multifractal formalism [15, 16] 
have provided insights into the structure of the mixing 
pattern and its decay rate. Of particular importance in 
some of these studies is the idea that for time-periodic 
flows the spatial mixing pattern becomes persistent, in 
the sense that it repeats itself in time but with a decreas- 
ing overall amplitude of fluctuations. Time-persistent 
spatial patterns have been observed in numerical sim- 
ulations [3, 6, 9], as well as in dye homogenization exper- 
iments in cellular flows [17-19], and have been related to 
the slowest decaying eigenmode of the advection-diffusion 
operator. The term strange eigenmode, originally coined 
by Pierrehumbert [3], is used to describe these patterns. 



The eigenmode amplitude decays exponentially with time 
at a rate determined by its associated eigenvalue, and an 
exponential decay of concentration variance has indeed 
been observed in various systems [3, 6, 8, 9, 12, 17-19]. 

However, such results were obtained either in ideal- 
ized systems [3-6, 8, 9, 20], or in cellular flows [17- 
19]. It has been suggested [21-24] that mixing might 
be slower in large-scale bounded flows, because of slow 
stretching dynamics in the vicinity of a no-slip wall. The 
specific form of the velocity field at a no-slip wall was 
first noticed and exploited by Chertkov and Lebedev [21]. 
The authors calculated concentration statistics by en- 
semble averaging over different realizations of a flow in 
a bounded domain with random time-dependence. They 
obtained a transient algebraic decay of the scalar variance 
attributed to the influence of the wall, followed by an 
asymptotic exponential phase. Shortly thereafter, exper- 
iments of elastic turbulence in a microchanncl [25] showed 
an anomalous scaling of mixing dynamics with the Pcclct 
number, which the authors related to the predictions of 
Chertkov and Lebedev [21]. Detailed numerical simula- 
tions of scalar advection by a short-correlated flow in a 
bounded domain were recently performed by Salman and 
Haynes [24], who characterized the scalar decay with a 
multi-stage scenario that includes a transient algebraic 
decay. All theoretical and numerical studies [21, 22, 24] 
assumed that, in bounded flows, scalar fluctuations are 
rapidly completely exhausted in the bulk because of effi- 
cient stretching therein, while scalar inhomogeneity sub- 
sists only in a decreasing pool at the boundary. 

In a previous letter [26], we have reported on the first 
experimental observation of "slow" algebraic mixing dy- 
namics imposed by a no-slip wall in a deterministic 2-D 
chaotic advection protocol. In the present paper, we ex- 
plore in more detail the successive stages of mixing of a 
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passive scalar in experiments and simulations of Stokes 
flow. For several chaotic advection protocols, a blob of 
dye is released in a closed vessel and homogenized. In 
contrast with all studies mentioned above, we focus here 
on the influence of the wall on the concentration field far 
from the boundary, which we show is contaminated by 
the algebraic dynamics near the wall. Our approach is 
based on a Lagrangian description of stretched filaments 
slowly fed from the wall into the bulk. A simplified one- 
dimensional model — a generalization of the baker's map 
— allows us to describe the various mechanisms at play 
and to reproduce the main features of the evolution of 
the concentration probability distribution. 

The paper is organized as follows. In Section II we 
review the main ingredients of chaotic mixing. We dis- 
cuss the successive stages of mixing and the associated 
length scales, which are then illustrated on the peda- 
gogical example of the well-studied hyperbolic baker's 
map. For this ideal system, we relate the structure of the 
strange eigenmodc to the unstable manifold of the least- 
unstable periodic point. Section III is the core of the 
paper. We first report on homogenization experiments 
conducted with a figure-eight protocol already described 
in [26], which we complement by numerical simulations 
of a viscous version of the blinking vortex flow [1], and 
a modified version of the baker's map with a parabolic 
point at the wall. In all cases, we observe anomalously 
slow mixing, that is an algebraic — rather than expo- 
nential — decay of concentration variance. We argue 
that this behavior is generic for two-dimensional mixers 
where the chaotic region extends to fixed no-slip walls. In 
such systems, poorly stretched fluid escapes the wall at a 
slow rate (controlled by no-slip hydrodynamics) through 
the unstable manifold of parabolic points on the wall. 
These poorly-mixed blobs contaminate the whole mixing 
pattern, up to the core of the domain where stretching 
is larger. We show that the modified baker's map de- 
scribes the experiments qualitatively and allows an ana- 
lytic derivation of the observed scalings for the concen- 
tration distribution. A discussion on very long times, 
general initial conditions and hydrodynamical optimiza- 
tion is finally presented in Section V. 



II. HOMOGENIZATION MECHANISMS 
A. Stages and length scales of mixing 

In this section, we describe briefly how the concentra- 
tion field of a passive scalar (e.g. dye) evolves from an 
initially segregated state towards homogeneity. For illus- 
trative purposes, we will consider the example shown in 
Figure 1 of a blob of dye of initial scale 4>iob smaller than 
the velocity field scale £ v , which is of the same order as 
the domain size L. Very viscous fluids typically support 
only laminar flows. Such flows may still lead to compli- 
cated, that is chaotic, Lagrangian trajectories [1, 2]. 

Three different stages of the mixing process are ap- 



parent in Figure 1. An initial blob (Fig. 1(a)) is de- 
formed by the stirring velocity field. At early times, the 
concentration pattern evolves as finer scales are created, 
yet the variance (the spatially-integrated squared fluctu- 
ations from the mean) is almost unchanged as the spatial 
scales are still too large for diffusion to be efficient (first 
stage, Fig. 1(b)). After several stretching and folding 
events, the width of a filament of dye stretched at a typ- 
ical rate A stabilizes at the so-called Batchelor length 

w B ■= \Jk/\ (1) 

with k the diffusivity, where the effects of compression 
and diffusion balance. Obviously, in a realistic flow, 
stretching is not constant, but the width of a filament 
quickly adapts to the local stretching rate A(x). The 
length scale wb is the smallest that can be observed inside 
the concentration pattern: an initial blob with a scale 
greater than wb is stretched and folded into many fila- 
ments that are compressed up to the diffusive scale Wb- 
During a second stage (Fig. 1(c)), after a strip has stabi- 
lized at the width wb, the amplitude of the concentration 
profile decreases according to the stretching experienced 
by the strip, to ensure conservation of dye [27, 28] . Differ- 
ent gray levels correspond to different stretching histories 
along the elongated image of the initial blob. Finally, 
since filaments are stretched but also folded, they are 
eventually pressed against each other, and their diffusive 
boundaries interpenetrate (third stage, Fig. 1(d)). Ulti- 
mately, homogenization takes place inside a "box" of size 
Wb through the averaging of many strips that have expe- 
rienced different stretching histories and have therefore 
different amplitudes [29]. 

So far, the most satisfying explanation for the decay 
of inhomogeneity in chaotic mixing has been the strange 
eigenmode theory, initially proposed by Picrrchumbert in 
a 1994 paper [3]. The strange eigenmode is the second 
slowest decaying eigenmode of the advection-diffusion op- 
erator (the first trivial mode corresponds to a nondecay- 
ing uniform concentration) . We consider in the following 
the case of periodic velocity fields, where the periodically- 
strobed strange eigenmode is an eigenvector of the time- 
independent Floquet operator. It decays at an exponen- 
tial rate fixed by the real part of the corresponding eigen- 
value. The projection of the initial concentration field on 
this eigenmode decays slower than the contributions from 
other eigenmodes, so that one expects the concentration 
field to converge rapidly to a permanent spatial pattern 
determined by the strange eigenmode, whose contrast de- 
cays exponentially. 

A simple physical motivation for the strange eigen- 
mode is as follows. The asymptotic concentration is gov- 
erned by filaments that are pressed against each other in 
a box of width w B - But these filaments have explored the 
whole domain, and hence may possess different stretch- 
ing histories [30]. The decay rate can thus depend on 
global properties of the flow [9, 12]. In particular, it is 
sensitive to spatial correlations, and cannot be expressed 
simply in terms of the stretching statistics. 
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FIG. 1: Successive stages of homogenization for a blob of 
dye stirred by the figure-eight protocol (see section III for 
details), (a) An initial blob is stretched by gradients in the 
velocity field. At early times, stretched filaments are still 
too broad for diffusion to be noticeable, and the concentra- 
tion variance is constant, (b) As time increases, filaments are 
stretched and folded repeatedly, while a strip of white fluid 
coming from the boundary is inserted periodically at the core 
of the mixing pattern. As a result of this mass injection, the 
filamentary pattern grows slowly towards the boundary with 
time, while filaments become thin enough (c) for diffusion to 
become effective and cause the strips of dye to become more 
gray. Different gray levels correspond to different stretching 
histories, (d) Later, different filaments start interpenetrat- 
ing, and the concentration field results from the averaging of 
concentration values coming from neighboring strips. 



Since the seminal paper of Picirchumbert, strange 
cigenmodes have been observed in many numerical stud- 
ies [4-7, 9] . The proposed evidence for strange eigenmode 
were (i) the onset of permanent spatial concentration pat- 
terns and (ii) an exponential decay for the concentration 
variance, whose rate depended only weakly on the dif- 
fusion. Recurrent spatial patterns have also been ob- 
served in experiments where a viscous fluid is stirred by 
an array of magnets [17-19], yet the concentration decay 
seemed slower than exponential. In the following section, 
we briefly show how a strange eigenmode arises in a one- 
dimensional baker's map, and relate the spatial structure 
of the eigenmode to the regions of lowest stretching. 



B. Tracing out the strange eigenmode in a 
uniformly hyperbolic model 

In this section, we describe how periodic points of a 
uniformly hyperbolic map affect the concentration field 
C(x,t) of a low-diffusivity scalar, insofar as they deter- 



mine the spatial structure of the observed strange eigen- 
mode. We show that the concentration pattern obtained 
from an initial blob after successive iterations of the map 
is determined by the least unstable periodic point of the 
map, and its multifractal unstable manifold. For peda- 
gogical reasons, we use one of the most-studied model 
of chaotic mixing, the inhomogeneous area-preserving 
baker's map [4, 15, 16, 31]. 

The area-preserving baker's map is defined on a two- 
dimensional square region by dividing the region in two 
strips, stretching, and re-stacking them. It has the prop- 
erty of mapping a y-independent distribution to another 
such distribution [4]. We thus take our initial blob to be 
a strip uniform in the y-direction, and the baker's map 
stretches and folds this strip to create more strips, leaving 
the concentration independent of y. We can thus focus on 
one-dimensional distributions that depend only on the x 
coordinate: they represent a 'cut' across a striated pat- 
tern of strips like the pattern in Fig. 1. Hence, we limit 
ourselves to a one-dimensional version of the baker's map 
which captures the essence of dynamics. 

The baker's map / reads 

f:x^f 1 (x)Uf 2 (x) (2a) 

where 

h(x)=jx; f 2 (x) = 1 - (1 - 7 )z (2b) 

and the union (U) symbol in (2a) means that / is one- 
to-two: every point x has two images given by fi(x) 
and f2{x). The parameter 7 satisfies < 7 < 1 and 
controls the homogeneity of stretching, with 7=1/2 be- 
ing the perfectly homogeneous case. The baker's map is 
represented in Fig. 2. 

Under the action of the baker's map, the concentration 
profile evolves as 

C(x,t+l) = C(.r 1 (x),t). (3) 

/ therefore transforms the concentration profile C(x, t) 
at time t into two images "compressed" by respective 
factors 7 and 1 — 7. 

First, we compute numerically the evolution of an ini- 
tial blob under the action of /. The initial concentration 

C(x,0) = C (x) = J 1 ' x «^ x ^ x « (4) 
y 1 ' uw |0, otherwise, v ; 

is a strip of constant concentration between x a and x^. 
Diffusion is mimicked by letting the concentration evolve 
diffusively during a unit time interval [4, 30]. During that 
interval, C evolves according to the heat equation with 
diffusivity k. We use periodic boundary conditions. 

Figure 2 shows the concentration profile for a typical 
simulation with 7 = 0.4, after 17 iterations of / alter- 
nated with diffusive steps. We see that the system is 
well mixed, insofar as fluctuations of C around its spa- 
tial mean (C) (which is conserved by the map) are very 
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FIG. 2: (a) One-dimensional baker's map /, with two fixed 
points at x = and x = 1/(2 — 7). For 7 < 1/2 the most 
stable fixed point is Xf = 1/(2 — 7). (b) Concentration profile 
obtained for an initial blob transformed by 17 iterations of / 



(7 = 0.4, 



10 °). Dominant "spikes" are located at Xf and 



(with decreasing amplitude) around its iterates (of decreasing 
stability). 




FIG. 3: Main stages in the evolution of the concentration 
profile, from left to right. An initial blob is stretched into 
many filaments by the map. Once filaments reach the diffu- 
sive scale wb, intermediate concentration levels appear. The 
concentration profile takes the form of a persistent pattern — 
the strange eigenmode — when all boxes of size wb contain 
an image of the initial unit interval. 



weak compared to the initial blob. Angle brackets denote 
a spatial average. A closer inspection reveals that fluctu- 
ations of C have more important values at some points, 
so that the concentration pattern has distinctive spikes. 
Remarkably, the spatial pattern visible in Fig. 2(b) is 
permanent, as can be seen on the two rightmost profiles 
in Fig. 3: further iteration of the map does not change 
the form of C, only its amplitude. This shows that the 
concentration field converges very rapidly to an cigcn- 
function of the advection-diffusion operator, dubbed the 
strange eigenmode. We measure an exponential decay of 
the concentration variance, consistent with convergence 
to an eigenmode. Strange eigenmodes in the baker's 



map, and their decay rate, have been studied in detail 
in Refs. [4, 5, 20]. Here we provide a simple and novel 
way to characterize the spatial structure of the strange 
eigenmode. More explicitly, we describe below how the 
strange eigenmode pattern traces out the unstable mani- 
fold of the periodic point of the map with weakest stretch- 
ing. This echoes the description of invariant sets in open 
flows in terms of the chaotic saddle [32]. 

As can be seen in Fig. 2, the map / has two period-1 
(fixed) points, one at x — and another at x = Xf := 
1/(2 — 7). Of course, since the 1-D baker's map is one-to- 
two, both fixed points map to an additional point, so that 
the fixed points have iterates other than themselves. For 
7 < 1/2 (7 = 0.4 in Fig. 2), the second fixed point is less 
unstable than the first one, as the compression factor at 
Xf is smaller. We notice in Fig. 2 that the highest spike in 
the concentration pattern is located at x = X{, whereas 
spikes of decreasing height are located at iterates of xf. 
fi(xf), /2(/i(xf)), . . .: the unstable manifold of Xf forms 
the backbone of the concentration pattern. 

Let us examine more closely the iteration of the con- 
centration pattern. After t iterations of /, the initial blob 
has been transformed into 2* strips compressed by factors 
7*7 7*~ 1 (1 — 7): • • • 7 (1 — 7)'- However, diffusion imposes 
that the width of an elementary strip saturates at the 
Batchclor width wb where diffusion balances stretching. 
Using the range of possible stretchings in the map, we 
have 
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<wb< 



l-(l-7) 5 



where the diffusivity k has been rescaled by the size of the 
domain L = 1 and the time-period T = 1. We therefore 
approximate wb — \/k/(1 — T 2 ), where 

log r = 7 log(7) + (1 - 7) log(l - 7) 

is the Lyapunov exponent. The concentration profile of 
Fig. 2 has a typical variation scale of wb- Under re- 
peated compression and diffusion steps, each elementary 
strip converges to a Gaussian peak of width wb (see the 
second picture from left in Fig. 3), centered on iterates of 
the initial blob centroid x c := (x a +Xb)/2. The amplitude 
of each Gaussian strip is proportional to A/wb to con- 
serve total concentration, where A is the multiplicative 
compression experienced by the strip. 

The strange eigenmode regime is reached when the ini- 
tial blob has been stretched enough so that its centroid 
has an iterate in every box of size wb - that is, each box 
contains at least one image of the initial domain. In this 
regime, the concentration C(x,t) measured at a point x 
results from the addition of slightly shifted strips, whose 
centers are all iterates of x c that fit into a "box" of size 
Wb centered on x. The random averaging of such strips 
has been proposed as the mechanism controlling the ho- 
mogenization rate [29, 33]. However, due to strong time- 
correlations of stretching, this averaging is not a random 
uncorrelated process here. Fluctuations of C(x) around 
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(C) are typically inversely proportional to the number of 
iterates in the box centered on x [29] . High spikes in the 
pattern therefore correspond to "boxes" with relatively 
few contributing iterates — i.e. images of the initial do- 
main that have experienced relatively low compression. 

Iterates of the initial profile located around Xf have 
experienced successive 1 — 7 factors during the last steps 
of the process while converging to the attracting periodic 
point Xf , since they have been transformed by the second 
branch / 2 during all recent iterations. On average, these 
iterates have experienced compressions smaller than the 
mean compression T*. This explains why the sharpest 
fluctuations are visible around Xf, and with decreasing 
amplitude around iterates of Xf of decreasing stability. 
For example, all iterates around /i(xf) have experienced 
a large number of successive 1 — 7 compressions, and an 
additional 7 compression during the latest iteration. In 
contrast, iterates around Xf have experienced only (1 — 7) 
compression steps in the last iterations, and fluctuations 
are greater. 

We conclude that regions of low stretching control the 
structure of the concentration pattern. This behavior 
has already been illustrated for the case of mixed phase 
space (with elliptical islands or weakly connected chaotic 
domains) [7, 34], but also holds for purely hyperbolic 
domains with uneven stretching. 

In addition, our description of the coherent structure 
of the strange eigenmode sheds a new light on why it 
is so difficult to predict the decay rate of the eigen- 
mode [5, 8, 10, 20]. Fluctuations decrease at a rate de- 
termined by the subtle interplay of peaks corresponding 
to the successive iterates of Xf, therefore spatial correla- 
tions of stretching histories play an important role in the 
decay rate, which cannot be related easily to the distri- 
bution of stretching in the map. We have provided here 
an example of concentration patterns dominated by the 
periodic structures with least stretching. 

We now turn to the experimental study of mixing by 
fully chaotic flows in bounded domains. In such mixers, 
the dominant periodic structures are parabolic points on 
a no-slip boundary, and we show that such points impose 
slower algebraic dynamics. 



III. PARABOLIC POINTS AT THE WALLS 

In this section, we report on dye homogenization ex- 
periments conducted in a closed vessel where a single 
rod stirs fluid with a figure-eight motion. In this physi- 
cal system, the phase portrait is not purely hyperbolic as 
it was in the baker's map: we describe how separatrices 
(parabolic points) appear on the wall as a consequence 
of no- slip hydrodynamics. We show that these regions 
of low stretching slow down mixing and contaminate the 
whole mixing pattern up to its core, far from the wall. 
These experimental results were briefly presented in [26]. 
Here, results from a numerical simulation of a counter- 
rotating viscous blinking vortex protocol [1, 35] are also 



presented. The dye pattern bears a strong resemblance 
to that of the figure-eight protocol, and we show that 
parabolic points at the walls are again responsible for 
algebraic decay. Inspired by the baker's map studied 
in Sec. II B, we introduce a simplified 1-D model that 
produces comparable algebraic mixing dynamics for this 
broad class of mixers. 



A. Algebraic decay in experiments and numerical 
simulations 

We first describe the essential features of the experi- 
mental set-up [26] . A cylindrical rod periodically driven 
on a figure-eight path gently stirs viscous sugar syrup in- 
side a closed vessel of inner diameter 20 cm (Fig. 4(a)). 
The fluid viscosity v = 5x 10~ 4 m 2 -s _1 together with rod 
diameter i — 16 mm and stirring velocity U = 2 cm • s _1 
yield a Reynolds number Re = Uljv ~ 0.6, consistent 
with a Stokes flow regime. A spot of low-diffusivity dye 
(Indian ink diluted in sugar syrup) is injected at the sur- 
face of the fluid (Fig. 1 (a)), and we follow the evolution 
of the dye concentration field during the mixing process 
(Fig. 1). The concentration field is measured through 
the transparent bottom of the vessel using a 12-bit CCD 
camera at resolution 2000 x 2000. This protocol is a 
good candidate for efficient mixing: we can observe on 
a Poincare section (Fig. 4(b)) — computed numerically 
for the corresponding Stokes flow [36] — a large chaotic 
region spanning the entire domain, including the vicinity 
of the wall. 

We also perform numerical simulations of "dye homog- 
enization" for a different stirring protocol. We consider 
a viscous version of the blinking vortex flow [1, 35], that 
is, two counter-rotating vortices alternatively switched 
on and off. Following Jana et al. [35], we study a re- 
alistic version of this protocol consisting of two large 
fixed rods placed on a diameter of a circular domain (see 
Fig. 4(c)). To mimic the blinking vortex, the two rods are 
rotated one after the other through angles 9 and —6, in 
a counter-rotating fashion. This stirring protocol resem- 
bles the figure-eight, as the counter-rotating movement 
of the vortices draws fluid from the boundary in some 
part of the domain (the radial velocity v± is positive), 
whereas it is pushed towards the boundary in the other 
part (v± < 0). The flow parameters are 9 = 270° (angu- 
lar displacement of one rod at each half period), r = 0.7 
(distance between the rods), ri nnor = 0.2 (radius of the 
rods). Length scale units are irrelevant here, and all dis- 
tances are scaled by the radius of the cylindrical vessel 
f outer = 1 . In the same way, all times are rescaled by the 
stirring period T, so that we can work only with dimen- 
sionless quantities in the following. A Poincare section 
shows (Fig. 4(d)) that the chaotic region spans the en- 
tire domain for this protocol as well. The evolution of 
a "blob of dye" is mimicked by computing the positions 
of 10 6 particles — initialized inside a small square in the 
center of the domain — during 75 periods. 
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FIG. 4: [Color online] Homogenization in closed flows for the 
experimental realization of the figure-eight protocol (a)-(b) 
and the numerical simulation of the contra-rotating blinking 
vortex (c)-(d). The heart-shaped mixing patterns are very 
similar: an upper cusp corresponds to a parabolic injection 
point on the boundary, while in the lower part of the pattern 
filaments are nicely packed in a parallel fashion. Although an 
annular unmixed region remains in the vicinity of the bound- 
ary, the partly-mixed pattern grows towards the boundary 
with time. This purely chaotic phase portrait is confirmed by 
the Poincare section in (b) and (d), where a single trajectory 
fills the entire domain. In both cases two parabolic points can 
be inferred from the cusps in the upper and lower parts of the 
boundary. They correspond to separation points along whose 
unstable manifold fluid gets injected into the bulk, and to the 
corresponding reattachment at the opposite boundary. The 
frame in (a) and (d) indicates where measurements are taken. 
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FIG. 5: [Color online] Statistical properties of the concentra- 
tion field measured in a central region (see frames in Fig. 4 
(a) and (c)). The concentration variance (black circles) is 
consistent with the evolution law of Eq. (21) (gray solid line) 
- which is close to a power law - both in (a) the figure-eight 
protocol and (b) the blinking vortex. The decay of variance is 
much slower than an exponential, which would be the signa- 
ture of a strange eigenmode. (c) Concentration PDFs after 13, 
17 and 31 stirring periods in the figure-eight protocol. Both 
sides of the peak can be fitted by power laws (C ma x — C)~ 2 
(see solid line fits, and inset). Inset: left (light-gray) tail of 
the peak, P(C) against |C m ax(i) — C|. Also note the persis- 
tence of a white peak at C = 0, which transforms into large 
shoulder for longer times. 



In the experiments, we measure the concentration field 
inside a large rectangle (see Fig. 4(a)) far from the wall. 
We plot the resulting variance and probability distribu- 
tions functions (PDF) of the concentration in Fig. 5. 
We observe a decay of the concentration variance a 2 (C) 
that is much closer to algebraic than exponential, (see 
Fig. 5(a)). In fact, as we will see in Sec. IV, the decay 
is well approximated by log (t\v /t) jt 2 . This behavior 
persists until the end of the experiment (35 periods) , by 
which time the variance has decayed by more than three 
orders of magnitude. We also note that PDFs of con- 
centration have a wide shape, with power-law tails on 
both sides of the maximum (see Fig. 5(c)). Moreover, 
the PDFs of concentration are highly asymmetrical. A 
persistent white peak at zero-concentration values slowly 
transforms into a large shoulder at weak concentration 
values. This implies that the light-gray wing of the peak, 



corresponding to concentration values smaller than the 
most probable value, is more important than the dark- 
gray wing on the other side of the peak. Finally, the most 
probable value shifts slowly with time towards lower val- 
ues. 

For the simulations, we compute a coarse-grained con- 
centration field in a large area far from the boundary 
(rectangle in Fig. 4(c)). The coarse-graining (0.01 here) 
scale plays the role of the diffusive cut-off scale wq. 
Again, we observe an algebraic evolution of the concen- 
tration variance (shown in Fig. 5(b)) that is well fitted 
by the same decay law as in experiments. 
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B. Hydrodynamics near the wall 

These experimental and numerical results are not con- 
sistent with the exponential evolution of a single cigen- 
mode of the advection-diffusion operator. In order to 
understand these scalings, we first consider the various 
mechanisms at play. We will show that the observed slow 
mixing arises from a subtle combination of hydrodynam- 
ics, and the nature of the phase portrait at the wall. 

As can be observed on the Poincare section in Fig. 4(b), 
the chaotic region of the figure-eight protocol spans the 
whole domain, and no transport barriers are visible. (El- 
liptical islands can appear inside both loops of the figure- 
eight for a smaller rod, but for a large enough rod we did 
not detect such islands.) In particular, trajectories ini- 
tialized close to the wall boundary also belong to the 
chaotic region. They eventually escape from this periph- 
eral region to visit the remainder of the phase space, but 
only after a long time, as trajectories stick to the no-slip 
wall. This escape process takes place along the white 
cusp of the heart-shaped mixing region, as can be seen 
in Fig. 4. This white cusp is bisected by the unstable 
manifold of a separation point at the wall (upper red dot 
in Fig. 4 (b) and (d)). The manifold divides trajecto- 
ries reinjected from the left and right of the separation 
point. Since stretching is very weak close to the wall, 
fluid drawn into the heart of the chaotic region from the 
wall is poorly mixed at the moment when it is injected, 
as opposed to fluid that has spent some time there al- 
ready. Because we inject the initial blob of dye far from 
the boundary, poorly-stretched fluid injected from the 
boundary into the heart of the chaotic region consists of 
zero-concentration white strips that are interweaved into 
the mixing pattern (see Fig. 1 (b), (c) and (d)). 

We can better characterize such white strips in terms 
of hydrodynamics near the no-slip wall. Consider the 
velocity field v near the vessel boundary. The wall can 
be treated as locally flat, and we define local coordinates 
£11 and xj_ that denote respectively the distance along 
and perpendicular to the wall. No-slip boundary condi- 
tions impose v\\ = for x± = (on the wall) and the 
corresponding first-order linear scaling for small x±, 

W|| = -A(x||) xj_ + 0(x 2 L ), near the wall. (5) 

Note that we are modeling the net velocity field, as ev- 
ident in the Poincare sections in Figs. 4(b) and 4(d), so 
we ignore the periodic time-dependence. Incompressibil- 
ity implies 



dv\\ dv± 

+ — = 0, 



dx 



dxj_ 



which combined with (5) yields 



U'ix^xl + Oixi). 



(6) 



(7) 



Now from the Poincare sections Figs. 4(b) and (d) we can 
see that the only trajectories that consistently approach 



the wall do so along a separatrix connected to the wall 
in the lower part of the vessel (the lower dot in each fig- 
ure). All other trajectories recirculate into the bulk. The 
separatrix corresponds to a flow re-attachment point on 
the boundary [35, 37, 38], which we refer to as parabolic 
points. (All points on the boundary are parabolic fixed 
points, but the important ones for us have separatrices 
emanating from them. We only mean those distinguished 
points when we refer to parabolic fixed points.) 

If we choose x\\ = to be the position of the lower 
separatrix, then the velocity on the separatrix is 



0, 



\ax\. 



A'(Q), (8) 



since A(0) = in order that the separatrix and the wall 
be on the same streamline. Note that the linear part of 
the flow around the re-attachment point vanishes, hence 
the name parabolic. The requirement that particles ap- 
proach the wall along the separatrix implies a > 0. Inte- 
grating x± = v± using Eq. (8), we find 



x±(t) 



x 



1 + atxo ' 



(9) 



where x is the initial x± coordinate of the particle. 
Equation (9) predicts that the distance d(t) between the 
wall and a particle on the lower separatrix shrinks as 



d(t) ~ 1/at, t»l/ax - 



(10) 



This scaling was already derived in [21], from the same di- 
mensional reasoning. The rate of approach along the sep- 
aratrix constrains the approach to the wall of the entire 
mixing pattern. We verified both in the experiment and 
in the simulation that d(t) is indeed well-approximated 
by a power-law scaling d(t) ~ 1/t. Note also that Eq. (10) 
implies that particles along the separatrix 'forget' their 
initial condition for long times. This can be seen in 
Figs. 4(a) and (c): material lines 'bunch up' against each 
other in the lower part of the domain faster than they 
approach the wall. 

To ensure mass conservation, a quantity of unmixed 
white fluid scaling as d(t) oc t~ 2 is injected periodically 
in the mixing pattern. As each newly injected white strip 
has approximately the same length (determined by the 
extent of the rod path), the width A(i) of a strip injected 
at time t must also scale as 



A(i) := \d(t)\ = 1/at 2 



(11) 



where time t has been rescaled by the period T = 1 . 

The origin of our slow scaling now emerges. Clearly, 
the mixing patte lotically stretched and folded by 

the rod at each half-cycle, in the same manner as in a 
baker's map. Yet the folds are not stacked directly onto 
each other but are interweaved with the most recently in- 
jected white strip. Since each new white strip has a large 
width that decreases only algebraically with time, the de- 
cay of concentration is slowed down by this injection of 
unmixed material. 
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The dominant mechanism for mixing can be summed 
up as follows: (i) chaotic stretching imposes that the typ- 
ical width of a filament of dye in the bulk (i.e. far from 
the wall) shrinks exponentially down to the diffusion or 
measurement scale; yet (ii) wide strips of unmixed fluid 
of width A(t) oc t~ 2 are periodically interweaved with 
these fine structures. Both protocols have in common a 
chaotic region that spans the entire domain, which im- 
poses the presence of parabolic separation points on the 
boundary [35, 37, 38]. In the next section, we generalize 
the baker's map model to include such a parabolic point 
at the boundary, and reproduce the dominant features 
observed experimentally and numerically. 




(a) 



C. A modified baker's map model 

We now wish to derive quantitative predictions to ex- 
plain the observed algebraic scaling for the concentration 
variance. In the same spirit as in Section II B, we sim- 
plify the two-dimensional problem by characterizing only 
one-dimensional concentration profiles C(x, t) perpendic- 
ular to the stretching direction along which dye filaments 
align. The effect of the mixer during a half-period boils 
down to the action of a one-dimensional map that trans- 
forms concentration profiles by interweaving an unmixed 
strip of fluid with two compressed images of the profile. 
The width of each decays in time as A(i) oc t~ 2 , owing 
to the parabolic point on the boundary (Section III A). 
We therefore mimic the behavior of our mixer with a 
one-dimensional map g, in the spirit of the baker's map. 

The map g is defined on [0,1] for simplicity. It evolves 
concentration profiles as in Eq. (3) and satisfies the fol- 
lowing: (i) it is a continuous one-to-two function, to ac- 
count for the stretching and folding process; (ii) the 'wall' 
at x = is a marginally unstable (i.e. parabolic) point of 
g , so that the correct dynamics are reproduced by ex- 
panding g~ 1 (x) ~ x + ax 2 + ■ ■ ■ (a > 0), for small x; (iii) 
because of mass conservation, at each x the local slopes 
of the two branches, gi and g 2 , of g add up to 1. 

Other details of g arc unessential for our discussion. 
As in Section III A, diffusion is mimicked by letting the 
concentration profile diffuse between successive iterations 
of the map (with no-flux boundary conditions). This 
model is a modified baker's map, with a parabolic point 
at x = 0, as opposed to the baker's map in Section III A, 
where the dynamics are purely hyperbolic. The expres- 
sion of gi close to x = assures that the distance between 
the origin and successive iterates of a point by g\ shrinks 
as d(t) ~ 1/at, which is the same as Eq. (10) obtained in 
the experiments with no-slip hydrodynamics. 

We numerically evolve concentration profiles for the 
specific choice for g, 



gi(x) = x — ax 2 + (7 — 1 + a)x , 
g 2 (x) = 1 - ax 2 + (7 - 1 + a)x 3 , 



(12) 



with 7 = 0.55 and a — 
map 51(1) = g 2 {l) = 7, 



0.9. We fix k = 10~ 7 . In our 
and as for the baker's map we 
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FIG. 6: Concentration profiles from (a) the figure-eight exper- 
iment after 13 stirring periods, and from the 1-D model after 
(b) 25 and (c) 50 iterations of the map. Note the presence of 
"white fluid" (C = 0) near x = in all cases. 



approximate by 7 the mean stretching realized by g\ , and 
by (1 — 7) the mean stretching realized by gi — although 
stretching is not constant along the two branches. We 
choose 7 ^ 0.5 for an uneven stretching in the bulk, as 
in the experiments where fluid particles that stay close 
to the rod for long times experience more stretching 
than particles left behind. Our initial condition is of the 
form (4), with x a = (1 — S)/2 and x^ = (1 + 5)/2, i.e., a 
strip of width 8 centered on x c = 1/2. 

Figure 6 shows concentration profiles obtained after 



9 



several iterations of the map. Strong similarities are ob- 
servable between concentration profiles obtained in the 
experiment (Fig. 6(a)) and in the map (Fig. 6(b)-(c)). 
In both cases a thin layer of "white fluid" (C = 0) is 
present near x = 0. Its width decreases as I /at due 
to the parabolic point on the boundary. In the exper- 
iment, the concentration pattern in the bulk (far from 
the wall) is characterized by sharp spikes at zero or low 
concentration values, whereas fluctuations are quite weak 
elsewhere. The sharp spikes correspond to white strips 
recently injected from the boundary into the bulk. For 
the map, the bulk pattern (far from x — 0) is clearly 
dominated by a set of thin spikes, which are recently in- 
jected white strips. These white strips are images of the 
boundary region at x = by g 2 , which are successively 
iterated by g\ or g 2 after their injection in the bulk. 

The suitability of our model is also strengthened by 
statistical properties of the concentration field, which 
closely resemble the experiment. Figure 7(a) shows the 
concentration variance for the map (measured in a central 
region) superimposed with experimental data: again we 
find an algebraic evolution. Moreover, there is a strong 
similarity between the concentration PDFs depicted in 
Fig. 7(b) and the experimental ones shown in Fig. 5(c). 
In particular, they both exhibit power-law tails. We will 
see in Section IV that our modified baker's map is simple 
enough for the concentration statistics to be calculated 
explicitly 




(b) 



IV. CONCENTRATION STATISTICS FOR THE 
MODIFIED BAKER'S MAP 

The simplicity of the model introduced in Section III C 
allows us to calculate the statistical properties of the 
concentration field. In this section, our interpretation 
is based mostly on our simplified model, although com- 
parisons with the experiment are also made. We first 
consider the simple case of an initially-uniform blob, 
for which we characterize the concentration pattern by 
counting iterates of injected white strips, since these 
dominate the concentration pattern. We will treat more 
general initial conditions in Section V. We focus here on 
a central region where the concentration profile is at least 
partly mixed — that is, away from x = 0, where x is the 
map coordinate. The concentration PDFs and variance 
presented above have for example been measured in the 
range x £ [0.2, 0.9] in the map. We have checked that the 
variance measured in the whole domain evolves trivially 
as as it is dominated by the remaining white pool at 
the wall. (The variance in the numerical simulations of 
Salman and Haynes is measured in the whole domain and 
displays a t^ 1 evolution during the algebraic phase [24].) 
Here we are interested in the more complex evolution in 
the bulk, where stretching is high and the pattern seems 
"well mixed" after a few periods. 

The modified baker's map of Section III C transforms 
an initial blob of dye of width sq into an increasing num- 



FIG. 7: Variance and PDFs of the concentration field mea- 
sured in a central region for a blob of dye transformed by the 
modified baker's map (12). (a) The concentration variance 
(circles) shows an evolution close to a power-law, comparable 
to the figure-eight experiment (square symbols). The solid 
line fit and the t~ 4 slope correspond respectively to the evo- 
lution dictated by Eqs. (21) and (28). (b) The concentration 
PDFs have wide power-law tails on both sides of the peak. 
The dark-gray tail corresponding to high concentration val- 
ues decays much faster than the light-gray one (weak con- 
centration values), since the latter is continuously fed by the 
injection of white fluid from the boundary. 



ber of strips with widths sq Ai • • • A t , resulting from differ- 
ent stretching histories inside the mixed region, where A t 
is the compression experienced at time t. White strips 
also experience this multiplicative compression starting 
from their injection time. Because of diffusion, a strip 
of dye or white fluid is only compressed down to the 
local diffus ive Batchel or scale w B , that we approximate 
by wb = \Jk/{1 — T 2 ), where 



io g r 



log 



dg-\x) 



dx 



(13) 



is the Lyapunov exponent and (• • • ) is the spatial mean 
taken over the region of measurement. For a moderate 
stretching inhomogencity in the bulk, we expect T to be 
close to 1/2. In experiments as in simulations, we probe 
the concentration field on a pixel size, or box size, which 
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is smaller than wb- 

As in the baker's map, different values of C correspond 
to a different combination of superimposed strips in a box 
of size wb- We characterize P(C) by considering the dif- 
ferent widths of reinjected white strips that one can find 
in a such a box. We will distinguish between three generic 
cases corresponding to a partition of the histogram P(C) 
in three different regions (see Fig. 5(c) and Fig. 7(b)): 
a white (W) peak at C — corresponding to recently- 
injected white strips that are still wider than u>b, and 
light-gray (LG) and dark-gray (DG) tails corresponding 
to respectively smaller and larger concentrations than the 
peak (mean) concentration. Once we have quantified the 
proportion of boxes contributing to these different values 
of C, the variance is readily obtained as 

o\C) = J \c-{C)) 2 P{C)AC = ol + al G + al G . (14) 

We treat each region of the histogram in turn in Sec- 
tions IV A-IV C, and combine the results in Section IV D. 

A. White pixels 



From the definition (15) of the injection time ri(t), 

t-n ~ (log at 2 + log w B )/ log r (19) 

for large i, therefore 

n w ~ (log at 2 + log w B )/(at 2 log T). (20) 

The fraction of white pixels nw is plotted versus time 
in Fig. 8(a). (There arc no free parameters.) For later 
times, we find an excellent agreement between the data 
and the expression (20) for nw- Note that during the 
first few iterations nw is constant: this corresponds to 
the initial phase when dye strips are still wider than wb 
and diffusion is ineffective (i.e. up to t such that 6 x T* = 
wb)- The concentration variance is also almost constant 
during this initial phase, so we discard it. We deduce 
the contribution of the white pixels to the concentration 
variance for t < rw, 

^ n ir<\2 „ X 2 v logat 2 + logw B . . 
a w = n w (C) *6 x . (21) 

Of course, <r w = for t > rw , since by then there are no 
purely- white strips left. 



Let us start with white (zero) concentration values 
that come from the stretched images of white strips in- 
jected before time t. White strips injected at an early 
time have been stretched and wiped out by diffusion, 
that is their width has become smaller than %. A 
white strip injected at time to has been compressed to 
a width A(£o)A to+ i • • • A t at time t. We neglect the 
spatial variation of A(x) in the bulk and approximate 
A to+ i ■ • • A t ~ r*~ to . The oldest white strips that can be 
observed have been injected at time t — r;(t), where r;(t) 
is defined by 



A(7i)r*- T1 = t«B. 



(15) 



Note that t — r ; (t) - that is the number of periods needed 
to compress an injected strip to wb - is a decreasing 
function of time. After a time t — rw defined by 



A(r w ) 



wb , 



(16) 



the injected white strip is smaller than wb and no white 
pixels can be observed. For t < rw, we can observe all 
white strips that are images of strips injected between 
Ti(t) and t, and the number of white pixels is proportional 
to 

t 

n w = ]T A(n) = d(t) - d(7i(t)), (17) 

ra=Ti(t)+l 

where A(t) = d(t) -d(t+l). 

We now use the expression for d(t) imposed by dynam- 
ics close to the parabolic point, d(t) ~ I /at, which yields 
for large t 



B. Light-gray tail 

We now focus on the distribution of light-gray values 
corresponding to white strips that have just been com- 
pressed below the cut-off scale wb- A white strip is first 
injected between images of the mixing pattern where fluc- 
tuations are lower (see Fig. 6). Fluctuations measured 
in a pixel are therefore mostly due to a recently-injected 
white strip that is superimposed onto a homogeneous dis- 
tribution. We approximate the measured value C as the 
average of the largest white strip with width ^w < w b, 
and mixed "gray" fluid whose concentration is close to 
the most probable concentration C g . A box containing a 
white strip of scale £w thus carries a concentration 



C = C g x (1-^w/tife), 



(22) 



and we can relate the concentration PDF P(C) to the 
distribution Q(£w) of widths of images of the injected 
white strips in the following way: 



P(C) = Q(£ w ) 



d£w 



dC 



w B 



q(M- 



(23) 



Q(£w) is easily retrieved from standard combinatorial 
arguments. A white strip injected at to is transformed 
into 2*~*° images with scales A(i )r t_t ° (once again we 
consider only the mean stretching T, which amounts to 
matching a given concentration to a unique injection 
time). In a "quasi-static" approximation, we neglect the 
algebraic dependence of £ w (and hence of C as well) on 
to in the factor A (to) compared to the exponential de- 
pendence in r*~ to . Therefore 



nw - (*- ri)(arit)" 



(18) 



Q(£w) = (WA(to)) 



iog(2)/io g (r) 



x(l/* w logr), (24) 
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FIG. 8: Simulations of dye homogenization by the modified 
baker's map: (a) Fraction of white pixels (where C = 0) »w 
as a function of iterate (full circles). At early times diffusion 
has not yet started to broaden strips of dye, and % remains 
approximately constant. When nw begins to decay, it closely 
follows Eq. (20) (solid line) obtained by counting the images 
of injected white strips that have not yet been compressed 
down to the diffusion scale Inset: the distance between 
the dye pattern and the wall (measured by the position of the 
first peak in Fig. 6) evolves as d(t) = 1/at (solid line), (b) 
Probability of a light-gray concentration level with a given 
distance to the peak, i.e. probability that |C g — dC — C\ < e 
(C g ~ 0.1, dC = 3 x 10~ 2 , e = 10" 3 , as a function of time. 
P(C g — C) agrees well with the t~ 2 evolution predicted by 
Eq. (26) with F ~ 1/2 (solid curve). 



resulting in 



P(C) = [A(t)] lo ^ lo ^ r ^ -^f- 



w B 1 - 



c 



I"K 2 

io g r" 



= P(t)[C 8 -C] 



(io g 2/i og r)-i 



(25) 



P(C) thus has a power-law tail in the light-gray levels 
whose exponent depends on the mean stretching F. We 
observe satisfactory agreement between this prediction 
and both experimental data and numerical 1-D simula- 
tions (see Figs. 5(c) and 7(b)). Indeed, for the tail in 
Fig. 5(c) and Fig. 7(b) we measure P(C) ex. (C - C g )~ a 
with a < 2, consistent with T < 1/2, a fairly homoge- 



neous stretching, as expected for 7 = 0.55. Also note 
that the amplitude of the light-gray tail decreases with 
time as a power law, 

P(t) cx [A(t)\ log 2 / u *C 1 / r > oc t~ 2(log 2/ Mi/r)). (26) 

We have plotted in Fig. 8(b) the probability of a con- 
centration value at fixed distance from the maximum. 
The observed evolution scales as a power-law t~ 2 , as ex- 
pected from our calculation. We deduce the contribution 
of light-gray pixels to the concentration variance, 



'LG 



= P(t) / (C e ~C) 2 - a(r) dC, 



(27) 



where a(T) = 1 — log 2/ log T, and C m - m is the smallest 
concentration observed (C min = for t < tw and C m - ln — 
C g (l — A(t)/u; B ) for t > r w ). For t < r w the integral 

is constant and <Tlg x ^ x On the other 
hand, for t > ryy, 



J2 

°LG 



P(t) 

2 + a(r) 



ft-Cy^ai-W). (28) 



For a (r) ~ 2 as we observed, the exponent in the above 
power law is about —4. 



C. Dark-gray tail 

Let us now turn to the dark-gray part of the PDF. In 
our case, high concentration values correspond to black 
strips of dye that have experienced little compression, so 
that they have not been grayed-out by averaging with 
many other strips. This time, it is not sufficient to con- 
sider only the mean stretching T to characterize such 
strips as we did in Sections IV A and IV B, since stretch- 
ing histories far from the mean are involved. Looking at 
the concentration profiles in Fig. 6, we observe that the 
highest concentration values come from the reinjection of 
black strips pushed to the pattern boundary where they 
have experienced lower stretching than inside the pat- 
tern core. Such a positive concentration fluctuation is 
then mixed with the remainder of the pattern as succes- 
sive images are compressed by a factor of order T, in the 
same way as injected white strips. Many images of the 
initial blob may have aggregated inside a box of size Wb- 
If the decay of this highest-concentration "cliff" is slower 
than r* — the decay of an injected fluctuation inside the 
bulk — we can apply the same method for computing 
the shape of the dark-gray tail as we did for white strips 
and the light-gray tail. 

In the spirit of Eq. (22) , we write 



C = C g X (1 - 4lack/w B ) + 4lack/WB, 



(29) 



and relate the width 4iack to the injection time to as 
in Section IV B. This leads again to a power-law de- 
pendence (C — C g )~ 2 , this time for the dark-gray tail. 
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This is in good agreement with the observed scalings for 
both experimental and numerical PDFs (see Figs. 5(c) 
and 7(b)). 

We now wish to estimate the time decay of the am- 
plitude of this tail. To do so, we evaluate the amplitude 
of the highest concentration fluctuation, located on the 
left boundary of the pattern (see Fig. 6(b) and (c)). This 
will give us the concentration value for which a number of 
boxes of order 1 contribute to the histogram, and hence 
provide an approximation of the amplitude of the tail. 
The contribution of the dark-gray tail is tiny compared 
to the light-gray one, since, after a few periods, only a 
few boxes of size wb on the border of the pattern have 
an amplitude significantly greater than the mean (see 
Figs. 6(b) and (c)), whereas the width of the remain- 
ing white pool is much larger. Moreover, we show below 
that this amplitude decays faster than the contribution 
of white strips. 

We have plotted in Fig. 9 the decaying amplitude of 
the largest fluctuation in the pattern, that is of the left- 
most box in the mixing pattern (Figs. 6(b) and (c)). The 
evolution during 200 periods reveals a first exponential 
decay, whose rate increases with the diffusivity, followed 
by a power-law phase with an exponent of about —3. 
A simple analysis explains this evolution. The ampli- 
tude of this fluctuation can be estimated as A/w-q, where 

A ~ rii=i 5i(3i(l)) i s the compression factor experienced 
after t periods at the boundary of the pattern (i.e. by 
the leftmost blob image). (This is true as long as the dis- 
tance between the pattern and the wall is less than the 
diffusion scale at the boundary. We will discuss this final 
phase in Section VA.) For early times dye strips do not 
yet feel the effect of the wall, and the stretching factor 
can be approximated by 7* (we use 7 instead of T for 
evaluating the compression by repeated iterations of g\), 
and we expect the decay to be exponential with a rate 
log 7. This behavior is indeed observed for large enough 
diffusivitics (Fig. 9(a)). For small diffusivitics, few strips 
of dye are homogenized before the boundary of the mix- 
ing pattern reaches the wall region, where the effect of 
the parabolic point dominates, so we do not observe the 
first exponential phase. 

For long times g\(l) ~ {ai)~\ Si (#1(1)) - 1-2/i. The 
compression A can be approximated by 

t 

A~ 7 "° f\ (1-2/i). (30) 

i=n a + l 

The two factors in (30) account for (i) the exponential 
compression by successive factors of order 7 inside the 
bulk, and (ii) a weaker compression by factors converging 
slowly to 1 as the boundary of the pattern approaches the 
wall and experiences a compression determined by the 
parabolic point at x — (Section IIIB). The product 
rii=no+i(l — V*) converges to a power law t~ 2 for long 
times. The observed exponent is greater; this might come 
from a crossover between an exponential phase and the 
predicted t~ 2 phase. 
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FIG. 9: Amplitude C m ax of the maximum positive concentra- 
tion fluctuation, in log-linear (a) and log-log (b) coordinates, 
for different diffusivities. This corresponds to the amplitude 
of the leftmost spike of the concentration pattern, which has 
experienced weaker stretching than images of the initial blob 
in the core of the pattern. An exponential decay followed by 
a power-law evolution are evident. A line of slope —2 corre- 
sponding to the asymptotic evolution expected from Eq. (30) 
has been drawn for comparison. 



^From the above analysis, we see that positive concen- 
tration fluctuations decay fast with time, compared to 
the contribution of the remaining white layer at x = 0, 
which shrinks much more slowly. The contribution of the 
dark-gray tail to the concentration variance is very small 
compared to the light-gray tail, therefore we neglect it in 
the following computation of the variance. 

It is important to note that the asymmetry of the con- 
centration PDF persists because our mixer "remembers" 
the initial condition — a small black spot and a big white 
pool — even after very long times. 

In our experiments, additional contributions to the 
dark-gray tail come from dye particles trapped during 
some time in folds of the pattern where stretching is weak 
(notice the dark folds in Fig. 1(a)). This 2-D effect is not 
present in our map. The dark-gray tail therefore consists 
of contributions from the border of the pattern, but also 
from these folds. Nevertheless, we have checked that this 
contribution is small compared to the light-gray tail, and 
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decays rapidly with time. 

D. Total concentration variance 

We finally sum all contributions from different parts 
of the PDF to obtain <r 2 (C), as in Eq. (14). From the 
above discussion, we distinguish two phases, t < r\y when 
the variance is dominated by the contribution of recently 
injected white strips that have not yet reached %, and 
t > tw when the most important fluctuations come from 
the mixing of white strips and gray fluid. 

In the experiment, the crossover time rw is estimated 
as 30 periods. However, 3-D effects inside the fluid pre- 
vented us from conducting experiments for more than 
35 periods. For this early regime, fitting cr 2 (C) with 
er^ cx (21ogi + log u>B)/t 2 (gray line on Fig. 5(a)) gives 
good results, except close to rw where the contribution of 
the light-gray tail starts to dominate. On the other hand, 
in numerical simulations we observe (Fig. 7(b)) both the 
(21og£ + log u>b)A 2 behavior (black line), which can be 
interpreted as in the experiment, and the t~ 4 decay after 
t\v (100 periods for the case studied) given by c 2 G . For 
long times, the observed power-law arises from the spe- 
cific way of incorporating white strips whose width scales 
as t~ 2 inside the mixing pattern. 

Having established the origin of the scalings observed 
for the concentration variance and PDFs in our exper- 
iments, we turn in the next section to the analysis of 
asymptotically long times, different initial conditions, op- 
timization, and offer some concluding remarks. 

V. DISCUSSION 

We have explained in the preceding sections the main 
features of our experiment. In this section we tie two 
loose ends: we examine the long-time behavior of the con- 
centration in Section V A, and look at the effect of initial 
conditions in Section VB. Both aspects are more easily 
investigated in our simple map than in experiments, and 
our discussion is supported by numerical simulations of 
the 1-D map. In Section V C we address the issue of opti- 
mization of the mixing device based on what we learned 
about the role of walls. Finally, we close the paper with 
concluding remarks in Section VD. 

A. Recovering an eigenmode for long times 

We now consider the asymptotic regime, when we can 
no longer approximate the injected variance by the con- 
tribution of a white strip of width I /at 2 . This is because 
the mixing pattern is close enough to the wall that dif- 
fusion blurs the white layer at the boundary. For such 
large times, the mixing pattern can be described as an 
inverted half-Gaussian centered on x = (see Fig. 10(a)) 
that decays with time as fluid is reinjected in the bulk. 



At this time, fluctuations are very small in the rest of 
the pattern, and they are only controlled by the ampli- 
tude of the half-Gaussian. In this final regime, where 
the concentration pattern keeps a self-similar form with 
time, the concentration profile has eventually converged 
to an eigenmode of the advection-diffusion operator. The 
width of the half-Gaussian w is determined by the point 
where stretching and diffusion balance, 

w ° = v^3w ■ (31) 

Thus, with g'i(wo) = 1 — 2aw , we obtain for a small 
diffusivity 

w = (ft/4a) 1/3 . (32) 

We note that for small diffusivities, wo - the Batchclor 
scale at x — - is much greater than the Batchelor scale 
in the bulk Wb oc k 1 / 2 . 

Once the concentration at x = starts rising, which 
occurs for d(t) = 1 /at <~ u> (i.e. 350 periods for a diffu- 
sivity k= 10~ 7 !), the stabilized half-Gaussian decays ex- 
ponentially at a rate — log(g' 1 (w a )) — 2aw , which scales 
as k 1 / 3 . We have verified this scaling in numerical sim- 
ulations of the model (see Fig. 10(b)). Note that this is 
one of the very few examples where one can predict an- 
alytically the decay of a eigenmode (another noteworthy 
situation is the torus map considered in Rcf. [8]). How- 
ever, this eigenmode regime is not relevant in practice, as 
we only observe it when fluctuations are completely neg- 
ligible in the bulk. Its structure is also quite trivial: it 
consists of the half-Gaussian at x = 0, and of very small 
spikes centered on the iterates of x — in the bulk. 

The convergence to the eigenmode can be interpreted 
as follows. Once the mixing pattern reaches the diffusive 
layer wq at the boundary, every box with size equal to 
the local diffusive scale contains an iterate of the initial 
blob of dye. A global decay of the concentration variance 
is therefore possible from then on. 

B. Other initial conditions 

For the sake of completeness, we have performed nu- 
merical simulations for other initial conditions than a 
blob of dye. In particular, we have simulated two dif- 
ferent situations, a cosine profile, corresponding to an 
initial condition 

C(x,t = 0) = 1 + cos(4tt2;), (33) 

and a random profile where we attribute to each pixel 
a random value between and 1. This rapidly vary- 
ing profile is quickly smoothed everywhere on the local 
Batchclor scale. The main difference between the two ini- 
tial conditions may be assessed as follows. For the cosine 
profile, the initial scale of variation for the scalar field is 
much greater than the Batchelor scale at the boundary 
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FIG. 10: (a) Structure of the eigenmode: an inverted 
half- Gaussian of width wo decays exponentially at a rate 
— \og(g' 1 (wo)) = 2awo. (b) Concentration variance measured 
in the whole unit interval [0, 1] for different diffusivities. (The 
thick solid, dashed, solid, dot-dashed and dotted lines corre- 
spond respectively to the following values of k: 5.4 x 1CP 9 , 
1.3 x 1(T 8 , 4.3 x 1(T 8 1.2 x 10~ 7 and 3 x 10" 7 .) For long times, 
the evolution of the variance is exponential, corresponding to 
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FIG. 11: [Color online] Homogenization for two different ini- 
tial conditions: (a) C(x,t = 0) = 1 + cos(47r:r), (b) random 
initial profile. The scale of variation of the initial profile is 
much greater than wo in the first case, whereas it is of or- 
der wo close to the boundary in the random case. Main axes: 
concentration profile after 30 iterations of the map. Note that 
all variance is contained in the leftmost image of the unit in- 
terval (always transformed by gi), and around the iterates of 
x = 1. Inset: zoom on the border region for periods 30 (dark 
blue) to 180 (red), represented every fifth period. 



wq, whereas for the random initial condition, the scalar 
field already varies on the smallest possible scales. 

In the first case, as for the blob of dye, the scale of 
variation of the profile close to the boundary is large, of 
order g\(l) ~ (see inset of Fig. 11(a)). This case 

is therefore analogous to the blob of dye case. After a 
short time, most important fluctuations are concentrated 
in the leftmost image of the initial unit interval, that was 
iterated only by g%. Other iterates have wandered in the 
bulk where stretching is much more efficient, so that all 
fluctuations have died out — except for newly-reinjected 
iterates. The history of newly-reinjected iterates can be 
coded as a sequence G o g\ o g 2 o g± where G stands for 
the last few iterations, which corresponds to the reinjcc- 
tion inside the bulk of fluctuations at the left boundary. 
Even the leftmost iterate feels the spatial heterogeneity 
of stretching (see inset in Fig. 11(a)), as fluctuations ini- 
tially close to x = 1 have been more compressed, and 
they have overlapped and averaged. After some time, 



the profile at the boundary (inset in Fig. 11(a))) has a 
value significantly different from the mean concentration 
only at one or two "oscillations" , which is exactly what 
we observe for the blob of dye case (see Fig. 6). As in 
the latter case, we observe a power-law decay for the vari- 
ance evolution, which can be accounted for by the same 
reasoning. 

For the random profile case, the concentration profile 
at the boundary (inset in Fig. 11(b)) is much less coher- 
ent over successive periods. Indeed, the scale of variation 
of the concentration profile saturates immediately at the 
local Batchelor scale in the whole domain. The strips 
reinjected in the bulk result from the averaging of many 
strips at the boundary, and their amplitude is more dif- 
ficult to predict. We measured a non-monotonic decay 
of variance inside the bulk in this case, as the averaging 
of strips close to the wall depends on the instantaneous 
height of many neighboring strips (Fig. 11(b)). Yet, the 
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FIG. 12: (a) Evolution of the concentration variance in a 
large central domain for three different versions of the figure- 
eight protocol, corresponding to r-o = 0.24, 0.3 and 0.38, ob- 
tained by advecting Lagrangian particles in Stokes-flow nu- 
merical simulations, (b) Evolution of the concentration vari- 
ance for the same protocols vs time rescaled by the param- 
eter y/a (roo = 0.24, 0.3 and 0.38 correspond respectively to 
a = 0.39,0.72 and 1.39). A satisfying collapse of all curves is 
observed. 



strange eigenmode regime is only reached once fluctua- 
tions have died out everywhere, except for the leftmost 
box of size Wq and around the iterates of x — 0, so that 
there is a long transient phase also in this case. 



However, many features are common to all initial con- 
ditions that we checked. The spatial organization of the 
bulk profile is dominated by the unstable manifold of 
the parabolic point at x = 0, where fluctuations per- 
sist longer. As stretching is lower close to the bound- 
ary, the reinjected fluctuations are similar over succes- 
sive periods. The same reasoning as for the blob of dye 
yields concentration PDFs with power-law tails of the 
form \C — (C)|~ 2 , which arc indeed observed in all cases. 



C. Hydrodynamical optimization 

We have argued that for a wide class of mixing proto- 
cols the decay of the concentration initially obeys a power 
law. For industrial devices, it is of primary importance 
to optimize the decay during this initial phase, that is, 
to tune the prefactor in the power law. Our analysis in 
Section IV shows that the prefactor is essentially deter- 
mined by the parameter a, which controls the evolution 
of the distance between the mixing pattern and the wall 
in Eq. (10). (The exponent of the power law also depends 
weakly on the mean stretching T.) 

We check the validity of the paramctrization of the 
variance by the rate a for the figure-eight protocol. As 
for the blinking vortex protocol, we follow a large num- 
ber of particles, in order to record the evolution of a 
coarse-grained concentration field. We perform numeri- 
cal simulations of a Stokes-flow version of the figure-eight 
protocol, for different values of the radius of the figure- 
eight loops r^. In all cases, the initial condition is a 
small square of size 0.1 x 0.1 located close to the rod and 
containing 2.25 x 10 6 particles. We calculate the vari- 
ance of the concentration on a large half-crown of outer 
radius 0.8. We show results for the evolution of the vari- 
ance in Fig. 12 (a). For = 0.24, 0.3 and 0.38, we ob- 
serve a power-law evolution with an exponent close to 4, 
but slightly greater for the smallest radius — 0.24. 
Numerical simulations do not permit the same spatial 
resolution as in experiments; we use a coarse-graining 
scale wb = 10~ 2 . Measuring the distance of the mixing 
pattern to the wall, we deduce the values of a and rw- 
Using this coarse spatial resolution and for each value of 
a, we compute tw ~ 5. The evolution of the variance in 
Fig. 12 is therefore in agreement with the —4 exponent of 
the variance determined in Sec. IV for the regime t > rw- 
The larger exponent for the smallest radius could be at- 
tributed to a weaker mean stretching T in this case, since 
the rod travels in a smaller region. 

In our model, the dependence of the variance on the 
details of the protocol comes from the factor A 2 (i) ~ 
(y/at)~ 4 in the regime t > r w . We have plotted in Fig. 12 
(b) the evolution of the variance against the rescaled time 
\fat. We observe a very good collapse of all data on the 
same curve, supporting the idea that the main ingredient 
in the evolution of the variance is the parameter a. 

Increasing the mixing speed therefore amounts to in- 
creasing the rate a at which a particle approaches the 
wall. This can be achieved in a number of ways, such as 
increasing the rod diameter or by "scraping the bowl" , 
that is taking the stirrers closer to the wall as in Fig. 12. 
Determining the value of a as a function of hydrodynamic 
parameters such as the rod diameter or the size of the 
rod's orbit is beyond the scope of this article. However, 
our results suggest that in comparing different mixing 
protocols the rate a gives a simple estimate of the vari- 
ance decay rate and can replace more challenging mea- 
surements, such as the concentration field itself. 
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D. Conclusions 

The results of this paper can be summarized as follows. 
Using an approach based on the Lagrangian description 
of fluid particles stretched into filaments, we have high- 
lighted the role of least-unstable periodic structures in 
mixing dynamics, first in the well-known baker's map, 
and then for a broad class of 2-D closed flows where the 
chaotic region extends to a no-slip wall. For the latter 
class of systems, we have proposed a generic scenario for 
wall-dominated mixing dynamics. No-slip hydrodynam- 
ics in the wall region force poorly-mixed fluid to be slowly 
reinjected in the bulk along the unstable manifold of a 
parabolic point. (Note that phase portraits with many 
injection points are also possible. This does not affect the 
validity of our arguments.) Mixing dynamics are then 
controlled by the slow stretching at the wall, which con- 
taminates the whole mixing pattern up to its core. We 
observed a slow algebraic decay of the concentration vari- 
ance in experiments and numerical simulations, which we 
justified analytically using a 1-D model of a baker's map 
with a parabolic point on the boundary. An exponential 
decay corresponding to an cigcnmode is recovered in the 
model once iterates of the initial blob of dye are present 
in all boxes of size of the local diffusive scale, that is 
for extremely long times at which the variance has been 
almost completely exhausted. 

We characterize a mixing experiment by the follow- 
ing parameters: the flow's mean compression factor T; 
the algebraic rate a at which a particle approaches the 
boundary; the Batchelor scale % obtained from the dif- 
fusivity k and compression T; and the width of the initial 
blob 5. The successive stages of the mixing process inside 
the bulk can be summarized as follows. 

• For t < log(wB/5) /logT, all filaments are larger 
than wb and the variance is constant: <r 2 (C) = a 2 . 

— 1/2 

• For log (wb/S) /logT < t < (owe) , fluctua- 



tions in the bulk start to decay as dye filaments 
are compressed below wb, and the variance is dom- 
inated by large unmixed strips recently injected 
from the near-wall region into the bulk: a 2 cx 
S 2 x (log at 2 + log u! B )/(alogr x t 2 ), as derived in 
Sec. IV A. 

• For (oibb) -1 ' 2 < t < (n/Aa) 1 / 3 , all reinjected strips 
are smaller than wb, yet their contribution still 
dominates the variance evolution a 2 cx l/(a 2 t 4 ). 
This scaling was derived in Sec IV B. 

• For t > (k/Ao) 1 ^, we arc in the cigcnmode regime 
(see Sec. VA) and a 2 cx exp(— fit), where fi cx 
(m 2 ) 1/3 . 

Our study has highlighted the role of periodic struc- 
tures with lowest stretching in the construction of a time- 
persistent mixing pattern, dominated by their unstable 
manifold. This applies to the least-unstable periodic 
point in the baker's map (thus a hyperbolic point) and to 
the wall parabolic point for the figure-eight case. The im- 
portance of elliptic region for limiting mixing dynamics 
has been emphasized in other studies [7, 34]. 

In 2-D flows, where Lagrangian dynamics are Hamil- 
tonian, the wall region can either belong to a chaotic 
region, or to an elliptical island. We have argued that 
algebraic mixing dynamics are obtained in the first case. 
An experimental study of mixing dynamics in the second 
case is in preparation and will be reported elsewhere. 
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